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ABSTRACT 

Standard detection of evoked brain activity in functional 
MRI (fMRI) relies on a fixed and known shape of the impulse 
response of the neurovascular coupling, namely the hemo¬ 
dynamic response function (HRF). To cope with this issue, 
the joint detection-estimation (IDE) framework has been pro¬ 
posed. This formalism enables to estimate a HRF per region 
but for doing so, it assumes a prior brain partition (or parcel- 
lation) regarding hemodynamic territories. This partition has 
to be accurate enough to recover accurate HRF shapes but 
has also to overcome the detection-estimation issue: the lack 
of hemodynamics information in the non-active positions. 
An hemodynamically-based parcellation method is proposed, 
consisting first of a feature extraction step, followed by a 
Gaussian Mixture-based parcellation, which considers the 
injection of the activation levels in the parcellation process, 
in order to overcome the detection-estimation issue and find 
the underlying hemodynamics. 

Index Terms — joint detection-estimation, hemodynam¬ 
ics, Gaussian mixtures, parcellation, brain 

1. INTRODUCTION 

Functional MRI (fMRI) is an imaging technique that indi¬ 
rectly measures neural activity through the Blood-oxygen- 
level-dependent (BOLD) signal Q, which captures the vari¬ 
ation in blood oxygenation arising from an external stimula¬ 
tion. This variation also allows the estimation of the under¬ 
lying dynamics, namely the characterization of the so-called 
hemodynamic response function (HRF). The hemodynamic 
characteristics are likely to spatially vary, but can be consid¬ 
ered constant up to a certain spatial extent. Hence, it makes 
sense to estimate a single HRF shape for any given area of 
the brain. To this end, parcel-based approaches that segment 
fMRI data into functionally homogeneous regions and per¬ 
form parcelwise fMRI data analysis provide an appealing 
framework Q . 

In |[3}|^, a joint detection-estimation (IDE) approach has 
been proposed for simultaneously localizing evoked brain 
activity and estimating HRF shapes at a parcel-level. This 
spatial scale allows one to make a spatial compromise be¬ 


tween hemodynamics reproducibility and signal aggregation, 
the latter operation enhancing the inherent low signal-to- 
noise ratio (SNR) of fMRI data. This aggregation has to be 
well spatially controlled because of the known fiuctuation of 
hemodynamics across brain regions |[^[7|. Hence, the IDE 
approach operates on a prior partitioning of the brain into 
functionally homogeneous parcels, where the hemodynamics 
is assumed constant. A robust parcellation is needed to ensure 
a good IDE performance. A few attempts have been proposed 


to cope with this issue ||2||4||8f 11 

but most of them are either 

too computationally demanding 
for hemodynamics variability |[2 

1 

1Q||11|| or do not account 
M. Here, our goal is to 


provide di fast hemodynamics parcellation procedure that can 
be used in daily applications prior to IDE inference. 

In this work, we propose a two-step approach consisting 
first of hemodynamics feature extraction, in which a general 
linear model (GLM) is used to discriminate hemodynamics 
information, followed by a parcellation of these features. The 
goal here is finding features that are able to catch most of the 
hemodynamic information, without the need of perfectly esti¬ 
mating the HRF function. Afterwards, an agglomerative clus¬ 
tering algorithm is used to segment the features. Our main 
contribution is the consideration of the detection-estimation 
effect within the parcellation step: there is a lack of hemody¬ 
namics information in the non- or slightly-activating voxels. 
The idea is then to enforce grouping these uncertain voxels 
with neighboring activating voxels. This is done through a 
spatial constraint in the agglomerative step of the parcellation 
procedure. Moreover, the uncertainty in a given voxel can be 
quantified by a statistics linked to its activation level, namely 
a p-value obtained in the GLM feature extraction step. This 
statistics is hence injected within the agglomeration criterion. 

This paper is organized as follows. First, we present the 
artificial data generation, the feature extraction step in the 
GLM framework, and the proposed parcellation algorithm, 
which relies on a model-based agglomerative spatial cluster¬ 
ing. In the results, the proposed parcellation is compared with 
a classical Ward parcellation refined with a spatial constraint. 
Quantitative results are assessed in a Monte Carlo experiment 
where the variability against several random data sets is as¬ 
sessed. Finally, we investigate the impact of the input parcel¬ 
lation on the HRF estimation provided by the IDE approach. 








Fig. 1. Artificial fMRI data sets. Top left: hemodynamic 
parcels. Bottom left: response levels. Right: HRFs associated 
with parcels. 


2. HEMODYNAMICALLY BASED PARCELLATION 

2.1. Generation of artificial fMRI data sets 

Let us define the set of all parcels as P = {Vi^ ...^Vr} 

where is the set of position indexes belonging to parcel 
7 and denotes the set of positions in parcel 7. Artificial 
BOLD data sets are generated using the following regional 
BOLD model, for a given voxel j G and a given parcel 7: 

M 

yj = + Plj + bj , ( 1 ) 

m=l 

where is the response amplitude at voxel j for a cer¬ 
tain condition m, = {x'^_^^^)n=i:N 4 =Q:D is the binary 
matrix encoding the stimulus for each condition m, = 
{hdAt,j)d= 0 :D is the HRF corresponding to parcel At be¬ 
ing the HRF sampling period, P the orthogonal function basis 
multiplied by the drift ij , and bj the noise vector. Note that 
Vj = [yj{ti), where tk = k TR and TR > At. 
Typical values are At = 0.6 and TR = 1.8 sec. 

In this work, we considered artificial data at low SNR, 
with one experimental condition represented with a 20 x 20 - 
voxel binary activation labels q = and levels of ac¬ 

tivation a = [ai...aj], with {aj\qj = 1 ) ^ A/’(1.8,0.25). 
We simulated a map of hemodynamics parcels, with different 
HRF shapes (duration 25 sec., TR = 1 sec.. At = 0.5 
sec.) in each parcel (see Fig. [^, by using the combina¬ 
tion of 3 Bezier’s curves, each being controlled by 4 points, 
to describe the curves until the peak, from the peak to the un¬ 
dershoot, and from the undershoot to the end, given specific 
peak and undershoot widths. We considered a Discrete Co¬ 
sine Transform for P, a drift ij ^ A/’(0, IIJ 4 ), and white 
Gaussian noise with variance e| = 1.5. 
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Fig. 2. Hemodynamics features extracted 
Vj G J, and activation levels aj Vj G 


more fiexible regularized FIR (RFIR) approaches such as (T4) 
are also theoretically able to recover any HRF shape accu¬ 
rately in contrast to canonical GLM, RFIR inference is pretty 
difficult when the SNR is too low since it proceeds voxelwise. 
Moreover, it is time consuming because it performs unsuper¬ 
vised estimation (cf do))- Here, we are more interested in a 
quick feature extraction step that also allows us to disentangle 
true activated voxels from non-activated ones. 

As regards canonical GLMs, our feature extraction step 
consists in fitting the following linear model: yj = X[3jThj, 
where yj is the BOLD signal, A' the design matrix, (3j the 
parameter estimates and bj the noise at voxel j. Let us de¬ 
note ^j^o, Pj^i and pj ^2 the parameters associated with the 
regressors in A', which derives from the convolution of the 
experimental paradigm with the canonical HRF h, its tem¬ 
poral h' and dispersion derivatives h", respectively. We as¬ 
sume that /3j^o rather contains information about the HRF 
magnitude, whereas I3j^i and I3j^2 provide information about 
the HRF delay and dispersion respectively, and, hence, are 
useful to differentiate hemodynamic territories. Maximum 
likelihood (ML) inference enables to get the parameter es¬ 
timates f3j in each voxel among which we only retain (j)j = 
as input features to the parcellation method. 

To quantify the activation level, we consider the p-value 
Pj^o associated with testing Hq : Pj^o = 0 and we use the 
notation aj = 1 — pj^o G (0,1) for these statistics in voxel j: 
The higher the aj value, the larger our confidence in the pres¬ 
ence of evoked activity in voxel j. Importantly, the statistics 
aj does not enter in the parcellation along with the previously 
defined features cj)j , it is rather used as weights in the agglo¬ 
meration criterion. 

Fig.shows the features extracted Pjp (Fig.j^a)) and 

^,,2 (Fig.gb)), which contain information about the hemo¬ 
dynamics territories, as we can see in the activated regions of 
I3j^i with different values in the different territories; and the 
activation levels (Fig.|^c)) aj, with values from 0 to 1. 


2.3. Informed Gaussian mixture based parcellation 


2.2. Feature extraction 

As regards hemodynamic feature extraction, several ap¬ 
proaches are available. Here, we only focus on GLM-based 
ones involving either canonical HRF and its derivative(s) GD 
We chose not to rely in Finite Impulse Response (FIR) mod¬ 
els (Tg since they may be to sensitive to noise. Although 


2.3.1. Agglomerative clustering algorithms 

The model based interpretation of agglomerative clustering 
algorithms GD makes the assumption that features have been 
generated by probability distributions that vary across parcels. 
In the context of model-based hard clustering, which aims at 
assigning classes to the input samples instead of weights, the 
















goal is to maximize the classification likelihood with respect 
to both © and 2 : given a set of features 

j 

= ( 2 ) 

i=i 

where 0 = is the set of parcel-specific model 

parameters and 2 : = denotes the set of parcel labels 

associated with each voxel. 

In an agglomerative approach, this function is approxi¬ 
mately optimized by successive merge operations, starting 
from an initial clustering guess or singleton clusters. Hence, 
at each step, i.e. when merging two clusters Vj and Vr of 
the current parcellation 2 ; into the parcel V^' = U Vr 
of the resulting parcellation z', the relative increase of the 
log-likelihood has to be maximized: 

/max£( 0 ; \ 

logA£(©;^, 3 ;'|#) = log — ————— (3) 

I max£( 0 ; 2 : 1 / 

A given merging step thus involves several likelihood maxi¬ 
mization over parameters 0 . 

2.3.2. Gaussian-mixture model 

To account for the activation level aj associated with each 
voxel j, we rely on a independent two-class Gaussian mix¬ 
ture in the agglomerative step. The rationale is that features 
(pj are distributed differently within a given parcel depending 
on the corresponding activation levels aj. Hence, the two- 
class mixture is expressed on every in parcel 7 as a way to 
model parameter differences related to activation levels: 

1 

i=0 

1 

i=0 

The latent variable qj encodes the activation state of voxel j 
and Pr(g'j = 1) = refiects the probability of activation. 
This latent variable can be directly linked to the activation 
statistics aj obtained at the feature extraction step: = 

~ ^ “ ^ 7,0 where denotes the set of voxels 
in parcel 7 and their cardinality. Then, straightforward 
calculations give the following ML estimators for the parcel- 
level mixture moments: 

J GJ-7 

^7’* ~ ^ ^ A7,i)(0ji 1^1,i) !(^) 

J 

where = (1 - i - {-Ifaj), A = and 

and define the empirical weighted mean and covariance 
of features in parcel 7. Once the parameters 0 have been 


( 7^=5 

a 

Fig. 3. Parcellations obtained by SW and IGMM methods, 
considering an averaged 100 iterations Monte Carlo experi¬ 
ment. From left to right: noise variance 0,1, 2 and 5, respec¬ 
tively. 



Fig. 4. Ml-based quantitative comparison of SW and IGMM 
parcellation results with the ground truth territories for noise 
variances ranging from 0.0 to 5.0. 

estimated, the two parcels V^ and Vr that are selected to be 
merged into V^> = V^ UVr are those which maximizes: 

1 

logA£(©;5:,3;'|#) = X log(X X.i)) 

jGj-y i = 0 

1 

jelr i = 0 

1 

jei^, i=o 

3. SIMULATION RESULTS 

3.1. Informed Gaussian Mixture based Parcellation 

In this section, the proposed Informed GMM (IGMM) 
method is compared with the Ward fT^ algorithm with con¬ 
nectivity constraints, that here we name Spatial Ward (SW). 
A Monte Carlo experiment is used, where the variability 
against several random data sets is assessed, to quantify the 
results of both parcellation methods. Fig. shows results 
averaged across 100 runs and for different noise variance 
levels. SW makes a big parcel for all the non-activated po¬ 
sitions, while IGMM overcomes this issue and partitions the 
positions independently of the activation level. 
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Fig. 5. JDE detection results. From top to bottom: parcel- 
lation used, averaged estimated NRLs over 100 JDE itera¬ 
tions and detection MSB = being = [ef...ej], 

^true ^ [af and = aj - From left to 

right: ground truth, SW and IGMM. 


For a quantitative evaluation, mutual information (MI) (Tt) 
was used to compare our parcellation results with the ground 
truth territories. Fig. shows the evolution of both parcel¬ 
lation techniques with respect to increasing noise variance. 
IGMM outperforms SW and has a decreasing mean value 
until MI ^ 0.4, whereas SW has an almost constant mean 
MI « 0.25. Note that both methods are sensitive to noise 
and have a high variance. 

3.2. JDE analysis 

Here, we study the impact of the prior parcellation as input 
knowledge to the JDE procedure. As input parcellations, we 
considered the hard clustering resulting from either the SW- 
or IGMM-based average parcels computed over the 100 indi¬ 
vidual results of our previous Monte Carlo experiment, i.e., 
the ones shown in Fig. for = 1.5. 

Fig. compares the averaged detection results over 100 
JDE iterations. First, we can see a slightly lower mean de¬ 
tection MSB for SW-based hard parcellation. Also, a lower 
activation level can be observed for voxels considered in the 
wrong IGMM-based parcel. Fig. shows the averaged es¬ 
timated HRF profiles over 100 JDE iterations whether it is 
based on SW or IGMM parcellation methods, compared with 
the ground-truth HRFs in Fig. Overall, HRFs are well re¬ 
covered by both SW-based or IGMM-based JDE analyses. 
In region 2 (cyan), IGMM-based better fits the ground-truth 
since it mixes less voxels with different hemodynamics than 
SW which includes all non-activating voxels. In region 3 
(yellow), SW seems to yield a MC-averaged HRF estimate 
slightly closer to the ground-truth than IGMM which pro¬ 
duces a parcel that also spans region 4 (red). However, the 
MC variability is higher in this region than in the others for 
both parcellation methods as shown by the error bars. 

Regarding to computation time, IGMM takes 130 times 
more than SW, but 1800 times less than CC-JDE | [TT| . 


(a) Region 1. (b) Region 2. 
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Fig. 6. JDE HRF estimation averaged results over 100 JDE 
iterations, for the different regions, labelled in Fig.[2 The er- 
rorbars correspond to the standard deviation of the estimation. 


4. DISCUSSION AND CONCLUSION 

A hierarchical parcellation method that takes into account the 
activation levels in the parcellation process has been develo¬ 
ped, so as to be consistent with JDE assumption that one par¬ 
cel contains active and inactive positions, and find the under¬ 
lying hemodynamic territories, independently of the activa¬ 
tion level. In terms of computation load, it is not a demand¬ 
ing method. However, hierarchical agglomerative algorithms 
need the selection in advance of the number of clusters, and 
this fact lead us to a definition problem. In the JDE frame¬ 
work, we want parcels big enough to be able to estimate ac¬ 
curately the HRF, although if parcels are too big, we will be 
losing information and averaging the hemodynamics. How¬ 
ever, if we can have a good estimate of hemodynamics, we 
are interested in having as much as possible parcels to better 
recover the territories’ singularities. 

We quantitatively validated that the proposed IGMM ap¬ 
proach enables a better recovery than a reference spatial Ward 
approach. Indeed, parcels obtained with IGMM are less in- 
fiuenced by highly activated positions and do not mix non¬ 
activation positions altogether. Finally, JDE results are quite 
comparable in terms of detection and HRF estimation whether 
the input parcellation comes from the SW method or the pro¬ 
posed IGMM approach. Still, the proposed approach yields 
more reliable parcellations than SW and may be more ade¬ 
quate to treat real data sets, which will be investigated in fu¬ 
ture works. 
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